function EEG = cvc_load_bdf(bdffile,fs,elpfile)
% CVC_LOAD_BDF
%
%  Synopsis
%  ========
%
%  EEG = cvc_load_bdf(bdffile,fs)
%  EEG = cvc_load_bdf(bdffile,fs,elpfile)
%
%  -- Author: Mads Dyrholm --
%     Center for Visual Cognition, University of Copenhagen.
%     2009 - May 2013
%
%  Purpose
%  =======
%
%  Read Biosemi BDF file and create EEGLAB EEG-structure 
%  with sample rate decimation. Zero-phase anti-aliasing
%  lowpass filtering is used.
%
%  Inputs
%  ======
%
%  bdffile - Filename of BDF file.
%              
%  fs - New sample rate. Must be integer fraction of 
%  BDF sample rate, or it will be rounded to nearest 
%  integer fraction (at CVC the BDF sample rate is 2048Hz).
% 
%  elpfile - Location file for the EEG channels. Note:
%  only the EEG channels should be included in this file.
%  Bdf channels that match this file by name will get type
%  EEG. Other channels will be imported with the label from the
%  bdf header, and with no channel type.
%
%  Outputs
%  =======
%
%  EEG - EEGLAB structure.

if (nargin<3) elpfile = []; end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       maxrecs = 50;        % if you get OUT OF MEMORY, reduce this number !
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

edf = openbdf(bdffile);
trigchan = edf.Head.NS;
breakflag = (edf.Head.NRec > maxrecs);
if breakflag==1,
  BLOCKRECORDS = maxrecs;
else
  BLOCKRECORDS = edf.Head.NRec;
end
ALLEEG = [];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% decimation filter design.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
edfFs = edf.Head.SampleRate(1);
r = edfFs/fs;
if r~=round(r)
  r = round(r); % decimation
  fs = edfFs/r;  
  fprintf('\n\nCVC:  Warning!!!: Output sample rate adjusted to %f\n\n',fs);
end
nfilt = 8;
rip = 0.05;
[b,a] = cheby1(nfilt, rip, .8/r);
u = zeros(edfFs/2+1,1); u(round(length(u)/2))=1;
h=filtfilt(b,a,u);
delay = round(length(u)/2);
memory = [];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~isempty(elpfile)
    chanidx = []; 
    eloc = readlocs(elpfile);
    eloc(end+1).labels = 'TMP';
    for i = 1:size(edf.Head.Label,1)-1
        chanidx(i) = length(eloc);
        for j = 1:length(eloc)-1
            if strcmp(eloc(j).labels,strtrim(edf.Head.Label(i,:)))
                chanidx(i) = j;
                continue
            end
        end
    end
    EEG_chans = eloc(chanidx);
    for j=1:length(EEG_chans)
        if strcmp(EEG_chans(j).labels,'TMP')
            EEG_chans(j).labels = strtrim(edf.Head.Label(j,:));
        else
            EEG_chans(j).type = 'EEG';
        end
    end
end

% allocate
srate = edf.Head.SampleRate(1);
nsamp = edf.Head.NRec * srate;
downsampleddata = zeros(edf.Head.NS-1,ceil(nsamp/r),'single');
raweventchannel = zeros(1,nsamp);

% running indexes
dnidx = 1;
downphase = 0;

% get busy
fprintf('CVC: Downsampling on the fly, with anti-aliasing filter and group delay compensation ');
lin = 0;
for recordblock = 1:BLOCKRECORDS:edf.Head.NRec
  pct = round(100*recordblock/edf.Head.NRec);
  if floor(pct/10)>lin
    lin = floor(pct/10);
    fprintf('%i%%',pct);
  else
    fprintf('.');
  end
  
  [DAT,signal] = readbdf(edf,recordblock:min(recordblock+BLOCKRECORDS-1,edf.Head.NRec));
  EEGDATA = DAT.Record(1:end-1,:);
  pnts = size(EEGDATA,2);
  
  raweventchannel((recordblock-1)*srate+1:(recordblock-1)*srate+pnts) = bitand(uint32(16383),uint32(mod((DAT.Record(end,:)),2^24)));
  
  % correct somewhat for huge (!) DC offsets
  if (recordblock==1) 
    chanoffset = mean(EEGDATA,2);
  end
  for ch=1:size(EEGDATA,1)
    EEGDATA(ch,:) = EEGDATA(ch,:) - chanoffset(ch);
  end
  
  if r>1
    % filter.
    [EEGDATA,memory] = filter(h,1,EEGDATA,memory,2);
  
    % downsample the rest
    tmp = downsample(EEGDATA',r,downphase)';
  
    % correct for filter delay
    if (recordblock==1)
      tmp = tmp(:,max(round(delay/r),1):end);
    end
  else
    tmp = EEGDATA;
  end
    
  % insert
  updidx = size(tmp,2);
  downsampleddata(:,dnidx:dnidx+updidx-1) = single(tmp);
  dnidx = dnidx + updidx;
  
  % prepare for next downsampling
  downphase = mod(pnts-downphase,r);
  if (downphase>0) downphase = r-downphase; end

end

fprintf('\nCVC: extracting event latencies and types from event channel\n');
latencies = find(diff(raweventchannel));
types = raweventchannel(latencies+1);

%for ev = 1:length(latencies)
%  EEG_events.event(ev).latency = latencies(ev);
%  EEG_events.event(ev).type = types(ev);
%  EEG_events.urevent(ev).latency = latencies(ev);
%  EEG_events.urevent(ev).type = types(ev);
%end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%fprintf('\nDecimating the event latencies . . !\n');
%for ev = 1:length(EEG_events.event) 
%  EEG_events.event(ev).latency = max(EEG_events.event(ev).latency / r, 1);
%  EEG_events.urevent(ev).latency = max(EEG_events.urevent(ev).latency / r, 1);
%end

EEG = eeg_emptyset;
EEG.data = downsampleddata;
EEG.srate = fs;
EEG.ref = '';
EEG = eeg_checkset(EEG);

% (*)
EEG = pop_importevent( EEG, 'event', [types',latencies'], 'timeunit', 1/srate , 'append', 'no', 'fields', { 'type', 'latency' },'optimalign','off','align',nan);

% Quote from EEGLAB tutorial...
% -------------------------------------------------------------
% "latency - is the latency of the event. 
%  The information is stored in data samples (pnts) relative to
%  the beginning of the continuous data matrix (EEG.data),
%  but is displayed in the GUI in milliseconds relative 
%  to the beginning of the epoch."
% -------------------------------------------------------------

% Hmm, .it is actually 1-based. An EEG.event(-).latency of 1 gives an event
% on the first sample. I have checked this with pop_eegplot, and also
% importevent adds one due to the data sample-based indexing (see importevent code).

% Lets do it manually...
% Import the events
%eventstruct = importevent( [types',latencies'], [], 2048, 'timeunit', 1/2048, 'fields' ,  {'type'    'latency'}, 'append', 'no', 'align', nan)
% NB! importevent adds one due to the data sample-based indexing.
% so now, for decimating the events we will subtract it, hence the -1 in the following
% then add the 1 again at the decimated sample rate.
%for ev=1:length(eventstruct)
%  eventstruct(ev).latency = (eventstruct(ev).latency-1)/r + 1;
%end
% Check. This actually fits with pop_importevent above (*).
%EEG.event = eventstruct;

if ~isempty(elpfile),  
    EEG.chanlocs = EEG_chans;
    
    EEG = eeg_checkset(EEG);
    
    % fit locations to sphere
    omitchans = setdiff(1:length(EEG.chanlocs),eeg_chantype(EEG.chanlocs,'EEG'));
    EEG.chanlocs = pop_chancenter(EEG.chanlocs, [], omitchans); % nEEGdigitized+1:size(EEG.data,1));
    EEG = eeg_checkset(EEG);
end
EEG.ref = 'none';
EEG = eeg_checkset(EEG);

fprintf('\n DONE. \n');
